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Abstract 

The reduced basis method is a powerful model reduction technique 
designed to speed up the computation of multiple numerical solutions 
of parameterized partial differential equations (PDEs). We consider a 
quantity of interest, which is a linear functional of the parameterized 
PDE solution. Compared to the original quantity of interest, the quantity 
of interest computed using the reduced model is tainted by a reduction 
error. We present a new, efficiently and explicitly computable bound for 
this error, and we show on different examples that this error bound is 
more precise than existing ones. We also present an application of our 
work to certified sensitivity analysis studies. 

Keywords: reduced basis method, surrogate model, reduced order modelling, re- 
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Introduction 

A large number of mathematical models are based on partial differential equa- 
tions (PDEs). These models require input data (e.g., the physical features of the 
considered system, the geometry of the domain, the external forces...) which 
enter in the PDE as parameters. In many applications (for instance, design 
optimization, data assimilation, or uncertainty quantification), one has to nu- 
merically compute the solution of a parametrized partial differential equation 
for a large number of values of the parameters. In such a case, it is generally 
interesting, in terms of computation time, to perform all possible parameter- 
independent computations in an offline phase, which is done only once, and to 
call an online phase for each required value of the parameter, during which the 
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information gathered in the offline phase can be used to speed-up the computa- 
tion of an approximate solution of the PDE, and, hence, to reduce the marginal 
(ie., per parameter) computation cost. 

The reduced basis method |9] is a way of specifying such offline and online 
phases, which has been successuUy applied to various well-known PDEs |3j|6l 
[51118). One should note that, in the reduced basis (RB) method, the online 
phase does not compute a solution which is strictly identical to the numerical 
PDE solution, but an approximation of it, obtained by projecting the original 
discretized equations onto a well-chosen basis. In the application cases given 
above, however, one is not interested in the solution by itself, but rather in 
a quantity of interest, or model output, which is a functional of this solution. 
Taking this functional into account when performing the model reduction leads 
to a so-called goal-oriented method. For instance, goal-oriented basis choice 
procedures have been tried with success in the context of dynamical systems 
in [5ll20] , where the basis is chosen so as to contain the modes that are relevant 
to accurately represent the output of interest, and in a general context in pQ, 
where the basis is chosen so as to minimize the overall output error. All those 
papers showed that using an adapted basis could lead to a great improvement of 
reduction error. This paper is about goal-oriented error estimation, that is, the 
description of a rigorous and computable error bound between the model output 
and the reduced one. Two different reduced model outputs can be considered: 
the first (which we call uncorrected reduced output) is simply the output func- 
tional evaluated at the reduced output. The second (called corrected reduced 
output), described in |9], is the same, up to a correction term obtained from 
the solution of an auxiliary (dual) problem. The rate of convergence of the 
corrected output is better than the uncorrected one but the computation of the 
correction involves the application of the RB method to the dual problem, and 
this has generally the drawback of doubling offline and online computational 
times. Regarding output error estimation, ^ provides an error bound for the 
difference between the corrected reduced output and the original output. In this 
paper, we propose two new goal-oriented error bounds: one for the uncorrected 
reduced output, and one for the corrected reduced output. We also show, in 
numerical examples, that our bound is more precise than the existing bound. 
This paper is organized as follows: in the first part, we describe our output error 
bounds and explain how to compute them; in the second part, we see how to 
apply our error bound to certified sensitivity analysis studies; finally, the third 
and fourth parts present numerical applications. 

1 Methodology 
1.1 Preliminaries 

Reference problem We begin by setting up the context of the reduced basis 
method for afflne-parameterized linear partial differential equations presented in 
|9j. Our reference problem is the following: given a parameter tuple ji^V C M'f , 
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find the solution (in a discretized functional space X of finite dimension) 

of: 

Ai^i)ui^l) = fi^i), (1) 

where is a /^-dependent invertible square matrix of dimension dim and 
f dz X; then compute the output: 

si^i) = (2) 

where s : X — > K is a linear form on X . 

Important example of problem: Problems such as ([ij usually appear as dis- 
cretizations of /i-parametrized linear partial differential equations (PDE); the 
X space is typically a finite element subspace (e.g., Lagrange finite elements), 
and A{fi) and / are given by Galerkin projection of the weak form of the PDE 
onto a suitable basis of X. The boundary conditions of the PDE are usually 
either encoded in X or in A{^), and the inner product used to perform the 
Galerkin projection is typically the inner product. 

Choice of the inner product: We suppose that X is endowed with the standard 
Euclidean inner product: {u,v) — u^v, with associated norm ||w|| ~ ^ {u, u). In 
the PDE context described above, the use of the standard Euclidean product is 
justified by the fact that the relevant functional inner product (typically. L"^) has 
already been used to write A(/i) and /(/^), and the discrete matricial problem 
can be considered using the Euclidean inner product. 

Affine decomposition hypothesis: We suppose that A(/i) and /(/i) admit the 
following so-called affine decomposition [9j: 

Q Q' 

g=l 9'=1 

where Q,Q' £ W , Qq : V ^ Rand jq, : V ^ R {ior q ^ 1, ... ,Q, q' = 1, ... ,Q') 
are smooth functions, Aq are square matrices of dimension dimX and fq' G X. 
This hypothesis is required by the reduced basis method. 

Reduced basis method The dimension of the finite element subspace dimX 
is generally fairly large, so that the numerical computation of u(/i) from the 
inversion of A{fj,) is expensive. The reduced basis aims at speeding up "many 
queries", that is, the computation of u{fi) for all parameters ^ € Vq where Vq is 
a finite but "large" subset of the parameter set V. 

Reduced problem: We consider a subspace X of X, and a matrix Z whose 
columns are the components of a basis of X in a basis of X. This basis of X 
is called the reduced basis in the sequel. We denote by u(/i) the components, in 
the reduced basis, of the solution of the Galerkin projection of ([T]) onto X, that 
is, the solution of: 

Z'AZuip) = Z*/(m) (4) 
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(where, for any matrix M, M* is the transpose of M). 

Choice of the reduced subspace: There are different techniques for choosing the 
reduced basis (the Z matrix). This paper does not focus on this topic, but we 
cite the POD method (Proper orthogonal decomposition) [16], and the Greedy 
method [5]. 

Offline- online decomposition: The many-query computation can then be 
spht into two parts: the first part (usually called the "offline phase"), which 
is done only once, begins by finding a reduced subspace, then the Q parameter- 
independent matrices: 

Ag^Z^AgZ, q=l,...,Q 

and the Q' vectors: 

fq'^Z'f,,, q'^l,...,Q' 

are computed and stored. In the second part (the "online phase"), we compute, 
for each value of the parameter /i: 

_ Q _ _ Q' _ 

A{f,)^Y.®,{fi)A„ /(a^)- ^7,(m)/,' (5) 

q=l q' = l 

and solve for u(/.t) satisfying: 

A{^i)u{^,) = J{^l). (6) 

The key point is that the operations in ([5]) and ([6]) are performed on vectors and 
matrices of size dimX, and that the complexity of these operations is totally 
independent from the dimension of the underlying "truth" subspace X . In many 
cases, the smoothness of the map /i i~> u(/i) allows to find (in a constructive 
way, ie., compute) X so that dimX <C diiaX while keeping ||u(/i) — Z?I(/z)|| 
small, hence enabling significant computational savings. 

Output approximation: The output s(/^) can also be approximated from u(/i) 
using an efficient offline-online procedure: let Z e X be so that: 

s{u) = (Z, u) Vit G X] 

in the offline phase we compute and store: 

1= ZH 

and in the online phase we take: 

s(^) = (I, u{[i)) 

as an approximation for s(/i). 
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Reduced-basis error bounds 

Bound on u: Under additional coercivity hypothesis on A, the reduced basis 
method ^ also provides an efficient offline-online procedure for computing e"(/u) 
so that the approximation can be certified. This bound is based on the dual norm 
of the residual: 

pi^,)^\\A{p)Zu{^l)^f{^,)\\' 

where ||-|| is a suitably chosen norm on X (not necessarily ||-||), and a stability 
constant bound, which can be written as: 

a{n) > inf \v^A{fi)v\ (7) 

tiex,||D||'=i 

when A is symmetric. The inequality sign is due to the fact that the exact 
infimum can be costly to evaluate in the online stage; usually a procedure such 
as the successive constraints method ^ is used in order to quickly find a lower 
bound. 

The bound e"(/i) reads: 

The online procedure for the computation of e(/i) is also of complexity indepen- 
dent of c\m\X. 

Lipschitz bound on s: This online error bound can in turn be used to provide 
a certification on the output: 

y^ier \s{fj,)-Jifi)\ <m^e^ (8) 

and this bound is clearly "optimal" amongst those depending on /i through e^{^i) 
only. We call this bound the "Lipschitz"' bound, and denote it by e^(/i). 

The aim of Section [1.21 is to bound |s(/i) — s^(/i)| by a quantity which is smaller 
than (/i) of ([8]) and can be computed using an efficient offline-online procedure 
which does not require computation of e"(/i), described in Section Fl.BI In Section 
11.41 we consider a better approximation of s(/i) (denoted by sdlJ-)) which also 
depends on the solution of the adjoint equation of ([l} projected on a suitably 
selected dual reduced basis, and we see how the proposed bound for |s(/i) — 3^(/i) | 
can be modified in order to bound |s(/^) — Sc(m)I- 

1.2 Probabilistic error bound 

In this section, we give the expression of our output error bound. We begin 
with some notation: let's denote the residual by r{fi): 

r(/i) = A{^I)Zu{^l) - /(/i) e x, 
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and the adjoint problem solution (which will naturally appear in the proof of 
Theorem ll.il) by 

Let, for any orthonormal basis $ = {0i, . . . , ^tv"} of X, any N G N* , and 
i = l,...,N, 

We take a partition {Vi, . . . jVk} of the parameter space V, that is: 

'P = Uf=i-Pfc and k^k' ^ Vk n T'fe' = 0. 
We set, for i = 1, . . . , and k = 1, . . . ,K: 

^n.^ ($) = min A (m) , (*) = max A (a*) , 

and: 

p. [f^,'i>) j ^™n^($)else, 

\ else, 
where fc(/i) is the only fc in {1, . . . , K\ so that [i e T'fc. We also set: 

N N 

1 = 1 z = l 

Ti(a., A^, $) = max (|Tr (M)M7^r(A*)l) • 
Finally, we suppose that is a random variable on V and set: 



^(u;(^),0,)(r(^),0j 



i>JV 



where we take, for convenience, 0^ = for all i > A/" (so that the sum above is 
in fact between iV and M) . 
We have the following theorem: 

Theorem 1.1. For any a g]0; 1[ and for any N E N* , we have: 
p(^\s{^^)-^^^)\>T,{^i,N,<^>) + '^^^i^^ <a. 

Proof: We begin by noticing that: 

A{ny^r{n) = Zu{n) - u{n) 
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so that: 

s(/i) - s{fi) = {l,Zu{n) - u{n)) = {l,A{fiy'^r{fi)) = (w(/i),r(/i)). 
We expand the residual in the $ basis: 



i>i 



Hence: 



^^I)-s{^i)^Y.(^,A{^i)-'c|,,){r{^^),(t>,)=Y,H^^):<l>^){r{^^),^^)■ (9) 



We clearly have that for any N eE*: 

N N 



N 



i=l 



and this implies: 



N 



<Ti(/.,7V,$) 



(10) 



So we have: 



p[\si^,)-^^i)\>T,{^i,N,^) + 



r2(iv,$) 



<P|^|s(M)-s(/i)l > 



N 



Y,{r{^^)A^){A{^x)-%4>,) 

i=l 
N 

Y,{r{^^),^^){A{^i)-%<i^^) 



> 



a 

T2iN,<P) 



by (nni) 



< p 



J2{r{^^),<l>^){Ai^i)-n,q^, 



i>N 



> 



a 



by dl]) 



< a thanks to Markov's inequality. 



□ 



Choice of $ The error bound given in Theorem 11.11 above is valid for any 
orthonormal basis $. For efficiency reasons, we would like to choose $ so that 
the parameter-independent part T2(A^, <&) is the smallest possible, for a fixed 
truncation index N E N* . 

To our knowledge, minimizing T2{N,^) over orthonormal bases of X is an 
optimization problem for which no efficient algorithm exists. However, we can 
minimize an upper bound of T2{N, $). 
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We define an auto-adjoint, positive operator G : X ^ X by: 

1 
2 



ex, Gcl>=-E^{{r{fi),^)r{fi) + {w{fi),^)w{fi)). (11) 



Let Ai > A2 > . . . Xj\f > be the eigenvalues of G. Let, for i e {1, 2 . . . ,Af}, 
(j)f be an unit eigenvector of G associated with the i*'' eigenvalue, and = 

We can state that: 
Theorem 1.2. 

i>N 

Proof. We have: 

\i>Af i>Af / i>N 



Using Theorem 1.1 of [19J, we get that the minimum of T2^^{N,^) is attained 
for $ = ^'^ , and that minimum is J2i>N '-' 

This theorem suggests to use <i> = ^'^ , so as to control T2{N, $). 



1.3 Monte-Carlo approximation of the error bound 

In this Subsection, we present an implementable offline/online procedure for the 
estimation of the upper bound for |s^(/i) — presented in Theorem ll.il 

Estimation of (pf We fix a truncation index N E N*, and we estimate 
{<pf}i=i,...,N by using a modification of the method of snapshots used in Proper 
Orthogonal Decomposition [16] . This estimation is performed during the offline 
phase. We begin by estimating the G operator by G, then we approximate (f>f 
by the appropriate eigenvectors of G. 

Estimation of G: We take a finite (large), subset of parameters S C V, 
randomly sampled from the distribution of the parameter, and we approximate 
the G operator by: 

^^^^ E (('^(a^). <^)Km) + {wit,),^)w{f,)) 
" ties 

In other words, G is a Monte-Carlo estimator of G. We take {(/)p},;=i....,Ar as 
the unit eigenvectors associated with the N largest eigenvalues of G. 

Computation of the eigenvalues of G: The operator G admits the following 
matrix representation: 

G = [WW^ + RR*) , 
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where W (resp. R) is the matrix whose columns are the components of w{fj,) 

(resp. in a basis of X, for fi G 'E.. These two matrices have #S columns 

and dimX lines, which means that the matrix above is dhnX x dimX. 

In general, we take t^S dimX, and so it is computationally advantageous to 

notice that if </> is an eigenvector of G associated with a nonzero eigenvalue A, 

then: 

so that ip G Im W + Im R =: V. Hence, if V is the matrix of an orthonormal 
basis of V, then there exists ip so that (p = Vip and we have: 



WW^cp + RR*4) = X4> 



1 



2#S 



{{WW* +RR*)) V 



ip = Xip. 



As a consequence, it is sufficient to find the dominant eigenvectors ^p, . . . , i/)^ 
of the matrix S = -^VHWW* + RR*)V (of size 2S), and to deduce 0f from 



ipf by the relation 



Vipf ■ Besides, by writing E as: 



S = — ^ {{V*W){WW) + {V*R){R*V)) , 

it is possible to compute and store S without storing nor computing any dense 
dimX X dimX matrix. 



Computation of Ti (/i, N, $) All the quantities intervening in Ti {fi, N, $) can 
be straightforwardly deduced from ^™"''"'*^ and {r{fj,),(l)f). 

Computation of and (3™'^'' {<i>) : For i = 1, . . . , TV and fc = 1, . . . , if, 

the reals and /3[^°^(3>) can be computed during the offline phase, as 

they are parameter-independent. Thanks to the availability of the gradient of 
Di{fi) with respect to a quasi-Newton optimization such as L-BFGS [21] can 
be used so as to compute these reals; they can also be approximated by a simple 
discrete minimization: 

mm A(/i,$), ^IX^W- max A(A^,<i>), 

or, alternatively, we can compute the following parameter- independent, offline- 
computable quantities: 

(/,.,0f),(A,C,,0f) (z = l,...,A^, j = g-l,...,g, q' = l,...,Q') 

where {Ci, . . . , C„} are the column vectors of Z; they form a basis of the reduced 
space X. 

Computation of {r{ii)^(pf) : Let a parameter /i G be given, and let 

'Ui(/i), . . . , Un{tJ') be the components of the reduced solution u{fi) in the reduced 
basis {Ci, . . . ,Cn}- 
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By using the relation: 

Q n Q' 

9=1 J = l 9' = 1 

the dot products between the residual and <j)f can be computed in the online 
phase, with a complexity of 0{nQ + Q') arithmetic operations, 0{Q) evaluations 
of Q functions and 0{Q') evaluations of 7 functions, which is independent of 
dimX. 



Estimation of T2{N, $) We approximate T2{N, $) by computing the follow- 
ing Monte-Carlo estimator: 

N 

i=l 

As this quantity is ^-independent, it can be computed once and for all during 
the offline phase. 

Final error bound By using Theorem 11.11 we get that for e{^,a,N,^) = 
Tiip, N, $) + T2{N, $)/a, we have: 

P{\s{^I)-^^l)\>e{^I,a,N,^))<a, 

so we may take, as estimated (computable) error bound with risk a, 

?(M,a,jV,$)^ri(A.,jV,ci>)+ ' ^ (12) 

a 

In the rest of the text, this computable error bound is designated as the error 
bound on the non-corrected output, by contrast to the bound of the next part. 

1.4 Bound on the corrected output 

The reduced output s{^) is a natural reduced output that approximates s(/i). It 
is possible to solve an auxiliary problem in order to compute an error correction 
that improves the order of convergence of the reduced output. As we will see, 
the error bound presented above can be easily modified so as to certify the 
corrected output. 

Output correction The idea of solving an adjoint problem in order to im- 
prove the order of convergence of a scalar output has originated in [TU]. We 
introduced the so-called adjoint problem, whose solution itd(/i) satisfies: 



^(Ai)*u<i(A*) I, 
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and the solution UdifJ-) of the reduced adjomt problem: 

z'aA{fIYZdM^^) = z'j, 

where is the selected matrix of the reduced basis for the adjoint basis. 
The corrected output is: 

ScifJ.) = S{^) - {Zud{fi),r{n)) 

One shoud note that the corrected output can be computed using an effi- 
cient offline-online procedure which requires two reduced basis solutions, hence 
roughly doubling (when Z^ has the same number of columns than Z) the of- 
fline and online computation times, except in the particular case where A is 
symmetric and I is proportional to /. 

Existing error bound on the corrected output In [9], it is shown that: 

\s{\x)- Sc{\i)\< =: ecc(/i) (13) 

where a(/i) is the stability constant bound defined at ([7]) (in the symmetric 
case) and rd(/i) is the dual residual: 

rd(^) = A*{n)ZdUd{n) - Hp). 

Hereafter, the existing bound Scc is called dual-based error bound. 

Proposed error bound It is clear that the work performed in the above 
sections can be reused so as to provide a probabilistic bound on |s(/i) — Sc{fJ-)\, 
by simply replacing iy(/Lt) by: 

Wcin) = w{n) - Zudin), (14) 

and hence giving a competitor (called error bound on the corrected output) to 
the dual-based error bound. 

1.5 Summary of the different bounds 

To sum up, we have presented four output computable error bounds. Two of 
them are bounds for the uncorrected ouput error: 

• the "Lipschitz" error bound ([5]); 

• the estimated error bound on the uncorrected output that we propose 
in (ini); 

and two of them for the corrected output error: 

• the existing dual-based error bound ecc, defined at (jl3p : 

• the estimated error bound on the corrected output Cc, that is (|12p amended 
with p^ . 
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2 Application to sensitivity analysis 

Our error estimation method is applied in sensitivity analysis, so as to quantify 
the error caused by the replacement of the original model output by the reduced 
basis output during the Monte-Carlo estimation of the Sobol indices. For the 
sake of self-completeness, we briefly present the aim and the computation of 
these indices, and we refer to [TS], [13] and [7] for details. 



2.1 Definition of the Sobol indices 

For i = 1, . . . ,p, the i*'' Sobol index of a function of p variables s(/xi , . . . , /Xp) is 
defined by: 

g ^ Var {E{s{ni,...,fip)\fii)) 
Var(s(^i, . . . ,/Xp)) 

the variances and conditional expectation being taken with respect to a postu- 
lated distribution of the (/ii, . . . , fip) input vector accounting for the uncertainty 
on the inputs' value. These indices are well defined as soon as s G L'^i'P) and 
that Var (s(/ii, . . . , fip)) ^ 0. When fn, . . . ,fip are (stochastically) independent, 
the i^^ Sobol index can be interpreted as the fraction of the variance of the 
output that is caused by the uncertainty on the i^^ parameter in. All the Sobol 
indices lie in [0; 1]; the closer to zero (resp., one) Si is, the less (resp., the more) 
importance /ij's uncertainty has on s's uncertainty. 



2.2 Estimation of the Sobol indices 

The conditional expectation and variances appearing in ()15ll are generally not 
amenable to analytic computations. In those cases, one can estimate 5*^ by using 
a Monte-Carlo estimate: from two random, independent samples of size M of 
the inputs' distribution, we compute 2Af appropriate evaluations {sj} and {s'j} 
of s, and estimate Si by: 

^ ii Y.j=i ■SjSj - [if Y.j=i Sj] f if I]j=i 4) 

Si = ^ 5 —■ (16) 

When Ad and/or the required time for the evaluation of the model output are 
large, it is computationally advantageous to replace s by its surrogate model 
s. By using ([T6| on s (hence with reduced model outputs {sj} and {s'j}), one 
estimates the Sobol indices of the surrogate model rather than those of the true 
model. We presented in [7], Sections 3.1 and 3.2, a method to quantify the error 
made in the Sobol index estimation when replacing the original model by the 
surrogate one. We defined two estimators S"™^ and Sf\^_^ relying on 
output error bound samples {ej} and {Cj}, and proved that: 

Theorem 2.1. If: 

Vj' = 1, . . . , M, \sj — Sjl < tj and \s'j — 3^- 1 < Cj, 
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the 



have: 



In our case, the output error bound e{^) of Theorem 11.11 does not satisfy the 
above hypothesis, but satisfies a weaker "probabilistic" statement. This is the 
object of the following Corollary: 

Corollary 2.2. If: 

Vj = 1, . . . , M, P{\sj - Sjl > ej) <a anrf Vj = 1, . . . , M, P (|4 - s^^l > e'^) < a, 



then we have: 

P [S. e [Sl^^^j,; ^M-a„./2]) > (1 - X (1 - 
Proof. We easily have that: 



xP(Vj, Is,-?,! <e(Ai)) 
> (1 - a,,) X (1 - a)2^^ 



3 Numerical results I: Diffusion equation 
3.1 Benchmark problem 

Our benchmark problem |13) is the following: given a parameter vector 

M= (Ml, M2,A'3) eT'^ [0.25,0.5] x [2,4] x [0.1,0.2], 
we consider the domain — fl{fJ,) below: 



□ 



N 



Ml 



n 
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Our continuous field variable Ue = Ue{fJ-) € satisfies: 

Aue = in r2 
Ue = on T D 



|if = -1 on Tj, 

H = on \ (Fat U To) 



where 

Xe^{veH\n) s.t. «|r„=0}, 
A denotes the Laplace operator, and ^ is the normal derivative with respect 

to dn. 

This continuous variable denotes the potential of a steady, incompressible flow 
moving in a tube whose profile is given by J7, with open ends on Fjv and To- 
The Neumann boundary condition on Fjv states that the fiuid enters by Fjv 
with unit speed, the condition on d^} \ (Fat U F^) states that the velocity field 
is tangential to the boundary of the tube; finally the Dirichlet condition on 
Tj} guarantees well-posedness, as the potential field is determinated up to a 
constant. 

The problem (|17|) is equivalent to the following variational formulation: find 
Ue = Uein) S Xe SO that: 



VWe • Vti = — / V, Vv ^ Xe- 

n Jfn 

This variational problem is well-posed, as the bilinear form (m, v) i— > Vu • Vw 
is coercive on Xe (see, for instance, [17 , lemma A. 14). 

The above variational problem is discretized using a finite triangulation T of 
and the associated P^{T) (see or [T^) finite element subspace: find u E X 
so that 



Vu • Vw = - / V \/v e X, 

O Jfn 

where X = {v £ P^{T) s.t. vlr^ = 0}. 
In our experiments, dimX = 525. 

The affine decomposition of the matrix of the bilinear form in the left-hand side 
of the above equation is obtained by using a piecewise affine mapping from fl{p) 
to a reference domain D, as explained in |11) . page 11. 
Our scalar output of interest is taken to be: 

s{^I) = / u{^I), 

and V is endowed with the uniform distribution. 



3.2 Results 

We now present the numerical results obtained using the different error bounds 
on the output of the model described above. We report our bounds on the non- 
corrected and corrected outputs, as well as the dual-based output bound. Note 
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that the stabihty constant a{iJ,) is taken as the exact inf; this clearly advantages 
the dual-based output bound. 

For the comparisons to be fair, one should compare the error bounds of same 
online cost. It is widely assumed that there exists a constant C so that this cost 
is C X 2(dimX)'^ for the dual-based method, and C{dimX)^ for our method, 
since dual-based method involves online inversion of two linear systems of size 
dimX, and one system of the same size for our method. Hence, the reported 
reduced basis sizes for the dual method are multiplied by a factor \/2. 
In all cases, the reduced bases are computed using POD with snapshot size 80. 
To compute G, we use a snapshot of size 200. We also took K = 1 (ie., a trivial 
partition of V). The truncation index N is taken equal to 20. 
In Figure [U we compare the different error bounds on the non-corrected, and 
corrected output. For instance, for the error bound on the non-corrected output, 
we plot: 

where 5 is a random subset of V with size 200 and a, N, <&) is defined at 
(|12p . Other error bound means are computed accordingly. 

We also computed the mean of the Lipschitz error bound . It is not reported 
here as it was way higher than dual-based output error bound. We observe that 
our new output error bound outperforms the dual-based error bound, at least 
for finite reduced bases sizes. Two reasons can be seen to this superiority: the 
difference in nature (probabilistic vs. sure) between the two bounds, and the fact 
that we make a crucial use of expansion Q instead of using a Cauchy-Schwarz 
(or duality norm) argument. The rate of convergence (slope of the curve) of the 
corrected output is better than the non-corrected one, and this fact is reported 
by the two error bounds. Also, the expectation of T2 was estimated at 10~^^, 
which allows to choose a low target risk and remain very competitive. 



3.3 Application to sensitivity analysis 

We estimate confidence intervals for the sensitivity indices of s(/i) by using the 
method described in ff] , together with Corollary 12. 2[ and the non-corrected 
output. 

We take M = 1000 as sample size, B = 500 as number of bootstrap replications, 
dimX = 10 as reduced basis size, a = 0.00001 as output error bound risk, and 
aas — 0.05 as Monte-Carlo risk. The level of the combined confidence interval 
^rao./2; •^M-a„./2] then (1 - aa.)(l " a)*^ > 0.93. 

The results are gathered in Table [TJ The spread between 5™ and S^^ accounts 
for the metamodel-induced error in the estimation of the Sobol indices. The 
remaining spread between S*™^ ^^'^ the impact of the sampling 

error (due to the replacement of the variances in the definition of the Sobol 
indices by their empirical estimators). We see that, in this case, the metamodel- 
induced error (certified by the use of our goal-oriented error bound) is very small 
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3 4 5 6 7 8 9 10 11 12 

Equivalent reduced-basis size 



Figure 1: Comparison of tiie mean error bound on the non-corrected output, tlie 
mean dual-based error bound (tec) and tlie mean error bound on the corrected output 
(for risk a = 0.0001). The "equivalent" reduced basis sizes are in abscissae. 
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Input parameter 




cm . qM 

'-'i I '-'i 






cm . qM 




Ml 


[0.530352;0.530933] 


[0.48132; 0.5791] 


M2 


[0.451537;0.452099] 


[0.397962;0.51139] 


M3 


[0.00300247;0.0036825] 


[-0.0575764;0.0729923] 



Table 1: Results of the application of Section [2] to the estimation of the Sobol indices 
of the output of our benchmark model. 



with regard to the sampling error. We also notice that the estimate for the Sobol 
index for fi^ is negative; this is not contradictory as it is the true value of the 
index that is in [0, 1]. For small indices, the estimate can be negative. 

4 Numerical results II: transport equation 

We now apply our error bound on a non-homogeneous linear transport equation. 
Compared to the previous example, the considered PDE is of a different kind 
(hyperbolic rather than elliptic). 

4.1 Benchmark problem 

In this problem, the continuous field Ue = Ue{x,t) is the solution of the linear 
transport equation: 

— [x,t)+ t^-Q^Kx, t) = sm{x} exp(-a:;) 

for all {x,t) e]0, l[x]0, 1[, satisfying the initial condition: 
Ueix,t^O)^x{l~x) Va;e[0,l], 
and boundary condition: 

ue{x^o,t) = \fte [0,1]. 

The parameter fi is chosen in V = [0.5,1] and V is endowed with the uniform 
measure. 

We now choose a spatial discretization step Aa; > and a time discretization 
step At > 0, and we introduce our discrete unknown u — (u")i=o,...,Ar^;n=o,...,Art 
where 

1 1 

Nx = and Nt = — . 

Ax' At 

We note here that the considered PDE is hyperbolic and time-dependent, and 
that we perform the reduction on the space-time unknown u, of dimension 
{Nx + 1) • {Nt + 1). This is different from reducing the space-discretized equation 
at each time step. 
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The u vector satisfies the discretized initial-boundary conditions: 

Vi, u" = (iA.T)(l - il\x) (18) 

Vn, M^f = (19) 
and the first-order upwind scheme implicit relation: 

Vi, n ^ + — — = sin(zAa;) exp(-iAa;). (20) 

Ai Ax 

Let's denote by i? = -S(m) (resp. y) the matrix (resp. the vector) so that 
([lH]),([in]) and (Uni) are equivalent to: 



Bu^y (21) 

that is: 



B^Bu = B^y, (22) 

so that equation is ^ with A(/x) = B^B and / = B'^y. 
The output of interest is: s(/x) = uj^' . In the following, we take A< = 0.02 
and Ax = 0.05. As in the previous example, the true stability constants are 
computed for the dual-based error bound. 



4.2 Results 

We took a very low risk level a — 0.0001, a snapshot size of 200, iV = 10 retained 
(jif vectors and K = \. The results (Figure [5]) show that, once again, the error 
bounds we propose in this paper outperforms the dual-based error bound. We 
can also see the estimated error bound is always greater than the true error. 
This validates the empirical estimation of T2 by a Monte-Carlo approach, and 
it shows that the majoration of the risk level by a (i.e., the Markov inequality 
used to prove Theorem I l.ip is in fact quite conservative. 



Conclusion 

We have presented a new explicitly computable output error bound for the 
reduced-basis method. We have shown, on two different practical examples, 
that this bound is clearly better than the naive Lipschitz bound and that, at 
the expense of a slight, controllable risk, the performances of this new bound 
are better than the ones of the existing dual-based output error bound. 
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Figure 2: Comparison between the mean error bound (corr. and non-corr. outputs), 
the mean dual-based error bound, and the true (corr. and non-corr.) errors, for 
different reduced basis sizes. 
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